###
# 10/30/19
#
# Analysis of the image ratings to validate the SDVPS
#
library(tidyr)
library(dplyr)
library(lme4)
library(strengejacke)
library(psych)
library(pander)
library(knitr)

load("death-con-image-rating-data.RData")

# Pull out the "how disturbing items
howDisturbingDf <- ratingData %>% select(c(sh.1_3, sh.2_3, sh.3_3, sh.4_3, sh.5_3, sh.6_3, oh.3180_3, oh.6831_3, oh.9250_3, oh.9254_3, oh.9430_3, oh.9921_3))
howDisturbingDfL <- gather(howDisturbingDf, image, disturbing)
howDisturbingDfL$imageType <- factor(c(rep("sh",6), rep("oh",6)))
disturbingByImage <- howDisturbingDfL %>% group_by(image) %>% summarize(mean=mean(disturbing), sd = sd(disturbing))
disturbingByImage$image <- factor(disturbingByImage$image, labels=c("OH_3180", "OH_6831", "OH_9250","OH_9254",  "OH_9430", "OH_9921", "SH_1", "SH_2", "SH_3", "SH_4", "SH_5", "SH_6"))
disturbingByType <- howDisturbingDfL %>% group_by(imageType) %>% summarize(mean=mean(disturbing), sd = sd(disturbing))
disturbingByType$imageType <- factor(disturbingByType$imageType, labels=c("Other-Harm", "Self-Harm"))

# Drop the how-disturbing questions and then gather from wide to long format
ratingData <- ratingData %>% select(-c(sh.1_3, sh.2_3, sh.3_3, sh.4_3, sh.5_3, sh.6_3, oh.3180_3, oh.6831_3, oh.9250_3, oh.9254_3, oh.9430_3, oh.9921_3))
ratingDataLong <- gather(ratingData, "image_question", "rating", -studyID)

# Add in columns for the type of image, the imageID, and whether the value is realted to self harm or other harm
ratingDataLong$imageType <- factor(c(rep("sh",71*2*6),rep("oh", 71*2*6)))
ratingDataLong$imageID <- factor(c(rep("sh1", 71*2), rep("sh2", 71*2), rep("sh3", 71*2), rep("sh4", 71*2), rep("sh5", 71*2), rep("sh6", 71*2),
                                   rep("oh3180", 71*2), rep("oh6831", 71*2), rep("oh9250", 71*2), rep("oh9254", 71*2), rep("oh9430", 71*2), rep("oh9921", 71*2)))
ratingDataLong$questionType <- factor(c(rep("relatedToSH", 71), rep("relatedToOH", 71)))
ratingDataLong <- select(ratingDataLong, -image_question)

# We need to "widen" it a bit but it is still in long format for the analysis
ratingDataLong <- spread(ratingDataLong, value=rating, key=-c(studyID))

# Calculate related ratings based on image type
avgSHOHRatings <- ratingDataLong %>% group_by(imageType) %>% summarize(avgSH=mean(relatedToSH), sdSH=sd(relatedToSH),avgOH=mean(relatedToOH), sdOH=sd(relatedToOH))
avgSHOHRatings$imageType <- factor(avgSHOHRatings$imageType, labels=c("Other-Harm", "Self-Harm"))
kable(avgSHOHRatings, digits=2, col.names = c("Type", "mean-SH", "sd-SH", "mean-OH", "sd-OH"))

# Pull out the sh/oh ratings for self-harm images, run t-test and compute cohen's d
shImageRatings <- ratingDataLong %>% filter(imageType == "sh")
tResultsSH <- t.test(shImageRatings$relatedToSH, shImageRatings$relatedToOH)
pander(tResultsSH)
cohensDSH <- t2d(tResultsSH$statistic, length(shImageRatings$relatedToOH)*2, n1=length(shImageRatings$relatedToOH), n2=length(shImageRatings$relatedToSH))
sprintf("Cohen's D for Self-Harm Images rated as related to self-harm vs other-harm = %1.2f",cohensDSH)

# Pull out the sh/oh ratings for other-harm images, run t-test and compute cohen's d
ohImageRatings <- ratingDataLong %>% filter(imageType == "oh")
tResultsOH <- t.test(ohImageRatings$relatedToOH, ohImageRatings$relatedToSH)
pander(tResultsOH)
cohensDOH <- t2d(tResultsOH$statistic, length(shImageRatings$relatedToOH)*2, n1=length(shImageRatings$relatedToOH), n2=length(shImageRatings$relatedToSH))
sprintf("Cohen's D for Other-Harm Images rated as related to Other-Harm vs Self-Harm = %1.2f",cohensDOH)


# Perform straight regression model without accounting for repeated measures 
ratingFixedModel <- lm(relatedToSH~relatedToOH+imageType , data=ratingDataLong)
summary(ratingFixedModel)
tab_model(ratingFixedModel)

  
# Perform multilevel model taking into account only the repeated measures of the individuals
ratingMLMModel1 <- lmer(relatedToSH~relatedToOH+imageType  + (1| studyID), data=ratingDataLong)
summary(ratingMLMModel1)
tab_model(ratingMLMModel1, show.se=TRUE, show.stat=TRUE)
plot_model(ratingMLMModel1, type="re")

# Perform multilevel model taking into account the level two variable of imageType
ratingMLMModel2 <- lmer(relatedToSH~relatedToOH+imageType  + (imageType| studyID), data=ratingDataLong)
summary(ratingMLMModel2)
plot(ratingMLMModel2)
tab_model(ratingMLMModel2)

plot_model(ratingMLMModel2, type="re")

# Perform multilevel model taking into account the level two variable of imageType (flipping IV/DV above)
ratingMLMModel3 <- lmer(relatedToOH~relatedToSH+imageType  + (imageType| studyID), data=ratingDataLong)
summary(ratingMLMModel3)
plot(ratingMLMModel3)
tab_model(ratingMLMModel3)

plot_model(ratingMLMModel3, type="re")

save.image("death-con-sdvps-analysis-data.RData")

